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1 ABSTRACT 


The paper presents a time-domain method for computation of sound radia- 
tion from aircraft engine sources to the far-held. The effects of nonunifornr 
how around the aircraft and scattering of sound by fuselage and wings are ac- 
counted for in the formulation. The approach is based on the discretization 
of the inviscid how equations through a collocation form of the Discontinu- 
ous Galerkin spectral element method. An isoparametric representation of 
the underlying geometry is used in order to take full advantage of the spec- 
tral accuracy of the method. Large-scale computations are made possible 
by a parallel implementation based on message passing. Results obtained 
for radiation from an axisymmetric nacelle alone are compared with those 
obtained when the same nacelle is installed in a generic configuration, with 
and without a wing. 
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2 INTRODUCTION 


One of the main sources of aircraft noise is the engine. Most of the measure- 
ments and modeling regarding engine noise have been done for an isolated 
engine [1, 2]. The far field engine noise of an aircraft is, however, influenced 
both by the flow around the wings and fuselage and by scattering of sound 
from aircraft surfaces. Ideally, one would like to be able to compute the radi- 
ated noise field that takes into account these effects. One of the motivations 
behind this is to explore the possibility of reducing the aircraft noise foot- 
print by taking advantage of engine and wing location and by manipulating 
the nonuniform mean flow around the fuselage. This paper presents a spec- 
tral method suitable for large-scale simulations of engine noise propagation 
in a realistic configuration and flow environment around an aircraft. 

Recent advances both in computer architecture and computing methods 
have already prompted researchers to address the fan tone noise problem 
using Computational Aeroacoustics (CAA) time-domain techniques [3, 4, 
5]. A treatment of the acoustic liner boundary condition, which is not 
obvious in this context since the liner impedance data is usually available in 
the frequency domain, was recently proposed by Ozyoruk and Long [6, 7]. 
Stanescu et al. [5] used a multi-domain spectral method for fan tone noise 
radiation and showed that accurate results can be obtained with a relatively 
small (between four and five, depending on the highest polynomial degree) 
number of points per wavelength. The spectral multi-domain method is 
geometrically flexible, as it only requires finite-element type grids that can 
be generated for usual aircraft configurations with commercially available 
packages. Although this fact makes practical three-dimensional problems 
tractable, the computer resources needed for frequencies typical of those 
encountered in fan noise are still very large. Such problems can only be 
solved if a proper parallelization strategy is used in order to make best use 
of the modern distributed memory computers available. 

The discontinuous Galerkin (DG) method for hyperbolic systems of con- 
servation laws has been studied by Cockburn and Shu [8] in detail. A 
spectral element implementation was proposed by Kopriva et al. [9], and 
several studies have been performed on its dissipation and dispersion prop- 
erties [10, 11, 12]. These studies demonstrate clearly its advantages in terms 
of accuracy when applied to wave propagation problems. The DG method is 
otherwise completely similar to the multi-domain spectral method [5] both 
in formulation and implementation, and has therefore been considered as 
the discretization technique for this work. Through a proper choice of the 
flux computation points, the method only requires communication between 
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elements that have common faces, and is therefore very well suited for par- 
allelization using message passing. In this setting, an initial ordering of the 
elements pertaining to each processor allows communication of data for el- 
ements on partition boundaries to be overlapped with computation of the 
fluxes for those elements inside each partition, thus minimizing the impact 
of communication on the overall performance. 

The present paper is organized as follows. The next section briefly de- 
scribes the DG method and its numerical implementation. Issues related 
to geometry representation, grid generation and parallelization strategy are 
all addressed. The numerical results section then presents computations of 
the sound field generated by a single spinning mode propagating around 
an axisymmetric nacelle, and compares the results with those obtained for 
propagation of the same mode when the nacelle is mounted in a generic 
fuselage-pylon-nacelle configuration, both with and without a wing. The 
conclusions section ends the paper. 


3 THE DG SPECTRAL ELEMENT METHOD 


3.1 GOVERNING EQUATIONS 


The three-dimensional, nonlinear Euler equations of gas dynamics in con- 
servation form, 


9Q \r^ 9 F ( i 

dt ^ dx d 


= 0 , 


(1) 


are considered here to model the noise propagation process. The components 
of the state and flux vector are: 
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where p is the fluid density, Vi,i = 1,2,3 the fluid velocity components in 
the three Cartesian directions, p the pressure, E the total internal energy 
per unit mass, and Sij the Kronecker delta symbol. 


4 



3.2 DG DISCRETIZATION 


The computational domain is divided into non-overlapping general hexahe- 
dral elements that can have curved boundaries. Each element is then individ- 
ually mapped onto the master element Qm = [ — 1, l] 3 using an isoparametric 
transformation. Under the mapping, Eq. (1) becomes 


8Q *dF d 
dt + t x di d 


= 0 


( 3 ) 


where the new variables are the transformed components of the state and 
flux vector 

Q = JQ, F d = jJ2p^F m , (4) 

and J is the Jacobian of the transformation 


J = det 


d(xi,...,x 3 ) \ 
<9(£i, • • • ,£ 3 ) / 


( 5 ) 


The computational space coordinates will be denoted here by both (£ 1 , £ 2 , £ 3 ) 
and (£, 77 , C) for convenience. 

Let the space of polynomials of degree N in £ e [—1, 1] be denoted by Pjy. 
A basis for this space can be constructed using the Lagrange interpolating 
polynomials 

* 

hi(0= n 7T. ITT’ = (6) 

SjJ 

through the N + 1 Gauss-Legendre [13] quadrature nodes £j, * = 0, 1, . . . , N. 
The three-dimensional solution Q is searched in a trial space obtained as the 
tensor product of one-dimensional degree N polynomials for each coordinate 
direction, = Pn(£) x Pn(v) x R/v( C)> hence it takes the form 

N N N 
i=0j=0k=0 
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where Qijk(t) denote pointwise values of Q at time t. In view of a Galerkin 
discretization, the residual is required to be orthogonal to the same trial 
space within each element. Thus, the DG spectral element method looks for 
a solution of the system 


(Qii < fiijk ) + (V^ • F, 4> ijk ) — 0 (8) 

to be satisfied for all i,j,k = 0,1,..., N. Here (•,•) represents the usual 
L 2 inner product, and 4>ijk = hi(£)hj(ri)hk(() are the basis functions of Pjy. 
Using the divergence theorem, the equation is recast as 

{Qt, 4>ijk) + [ ■ ndS = (F, Vtfijk) (9) 

J dflM 

In order to solve these equations efficiently, the integrals are computed using 
Gauss quadrature 



N 

f(0d£ = /(Cm)^m, 

m = 0 


V/ e P2N+1 


(10) 


w m being the Gauss weights. This replacement is exact for elements with 
straight edges, in which case J is a polynomial of degree one. In the general 
case of elements with curved boundaries, a discretization error is incurred 
in the use of the numerical quadrature rule Eq. (10). 

Expanding the boundary integral and performing some algebraic manip- 
ulation, the final discrete form of the equations governing each variable at 
the Legendre-Gauss points can be found to be 


dQijk 

dt 


D* + D^ + D< 


F 


( 11 ) 


where the right-hand side is a sum of discrete differential operators acting 
on the flux values of an element, which include values on the element faces. 
Defining the contribution due to the volume term on the right-hand side of 
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Eq. (9) as 


N 

SF= Y;K j \Kihm)N 

771=0 


(12) 


N 

where the discrete inner product (u,v)n = ^ ~^u n v n w n was introduced, the 

71=0 

differential operator for £ for example can be written 


D^F 


1 


Wi 


F *0-,Vj,Ck)hi(l) 


F*(-l, Vj ,(k)hi(-l)-dSF 


(13) 


and the remaining two follow by obvious permutations. In this equation the 
notation F* denotes a common face flux as explained below. The position of 
the Gauss-Legendre points inside the element, as well as the position of the 
corresponding flux points on the element faces where F* must be evaluated, 
is shown in Figure 1, drawn for one ( = const plane and N = 2. 

It is important to note that the operation in Eq. (12) is actually a one- 
dimensional operation that is expressed in the computer code as the multi- 
plication of a matrix of entries di m = (h[, h m )N with a vector of size N + 1 
containing the values of the F\ flux component at the Legendre-Gauss points 
on a line i = const inside the element. These flux values can be computed 
directly from the values of the state vector which are readily available. The 
values of the state vector at the face points of each element, i.e. Q(l,rjjXk) 
and Q(—l,r]jXk) f° r the £ = const faces, are obtained by evaluating the 
interpolating polynomial based on the Gauss-Legendre points at £ = ±1, 
again a one-dinrensional (implemented as a vector dot product) operation. 
However, these values can not be used directly in computing the pointwise 
flux values on the element faces, as they are not uniquely defined. Indeed, 
for each face point, two different values are obtained for the value of the 
state vector Q, one from each of the two neighboring elements. In order 
to define a unique flux, a Riemann problem [14] is solved at each flux face 
point with the two different states as initial data, just as in the finite volume 
method. Letting the superscripts L and R denote the two (left and right) 
neighboring elements of a face, the application of the Riemann solver can 
be written 

*d(-) = F2(Q R (-),Q L (-)) (14) 

and completes the evaluation of the right-hand side of Eq. (11). Upon this, 
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it becomes on ordinary differential equation which is integrated in time using 
a low-storage Runge-Kutta scheme optimized for wave propagation [15] . 

The computer code thus obtained can be used for computing both a 
steady mean flow, when appropriate inflow/outflow boundary conditions are 
specified at the fan face and the far field, and the propagation of acoustic 
waves superposed on such a mean flow previously obtained. For acous- 
tic computations, once the mean flow is known, the boundary condition 
on acoustic source surfaces is allowed to become a function of time (cor- 
responding to a source term in the governing equations), and the acoustic 
pressure is obtained as the difference between the computed value of p and 
the mean flow pressure. Radiation boundaries are treated by a damping 
layer approach, similar to Stanescu [5]. The damping layer is activated only 
for acoustic computations. This precludes the need for a time-dependent 
radiation boundary condition, i.e. the boundary condition at the exterior 
limit of the damping layer is the same inflow/outflow characteristic bound- 
ary condition that can be used to compute the steady mean flow field. 


3.3 ISOPARAMETRIC ELEMENTS 

To obtain an accurate representation of the underlying geometry, the ICEM- 
CFD Hexa [16] commercial package has been used for hexagonal grid gen- 
eration around the configurations of interest. Once an unstructured grid of 
hexahedra was obtained, a new feature of this package allows for the gen- 
eration of points along each of the edges of the hexahedral mesh according 
to a prescribed distribution, which can be either a Legendre-Lobatto or a 
Chebyshev-Lobatto distribution in arclength, for a specified polynomial de- 
gree N. The resulting points are written to a data file which can be read in 
a computer code and used to define the geometry with the same accuracy 
as the solution. All the necessary point coordinates can then be computed 
by interpolation based on the spectral interpolants along the edges (to ob- 
tain coordinates at the Gauss points from the Lobatto points on the edges) 
followed by three-dimensional transfinite interpolation [17] on the faces and 
inside the elements. The intermediary use of the Lobatto points is found 
necessary because they include the eight corners of the hexahedra (and hence 
the end points for each edge of the grid) which can thus be guaranteed to 
have a unique position. This would not be the case if their coordinates were 
created by interpolation from data at Gauss points (they would be different 
due to interpolation and truncation errors). 



3.4 PARALLEL IMPLEMENTATION 

The DG spectral element method is particularly well suited for a parallel 
implementation through message passing as the computation is naturally 
structured by elements, with most of the work done inside the element i.e. 
evaluation of derivatives in Eq. (11). Here a brief description of the paral- 
lel implementation is given. The program starts with a preprocessing step 
wherein the underlying hexahedral mesh is decomposed in as many parti- 
tions as the number of processors using the multilevel graph partitioning 
techniques in METIS [18]. Each processor creates the connectivity for its 
grid partition, which includes a list of elements (’’cut elements” henceforth) 
having one or more ’’cut faces” in common with elements in neighboring 
partitions. 

One stage of the Runge-Kutta time integration scheme for the parallel 
DG algorithm can then be described as follows: 

• Compute state vector values at the face flux points through interpo- 
lation for all cut elements. 

• Post non-blocking send/receive for the state vector values on the cut 
faces. Each processor sends the values interpolated from the cut ele- 
ments it owns to the corresponding neighboring partition. 

• Compute state vector values at the face flux points through interpo- 
lation in the regular elements inside each partition. 

• Compute unique fluxes through the Riemann solver for regular faces, 
Eq. (14). 

• Test for completed send/receive. 

• Compute unique face fluxes through the Riemann solver for cut faces, 
Eq. (14). 

• Compute residual and update the state vector according to Eq. (11). 

The parallel performance of this algorithm has been tested on several 
architectures, for various grid sizes and polynomial degrees N. Results ob- 
tained for two grid sizes are presented here. The first one is a small grid 
of 884 elements with polynomial degree N = 12, which has been run on an 
IBM SP3 machine with 375-MHz Power3 processors, each with 512MB of 
memory. The second one is a medium-size (17370 elements) grid used for 
the fuselage-nacelle configuration, run with N = 4 on an SGI 0rigin2000 
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machine with 250MHz R10000 processors. Both cases could be run on a 
single processor, and Figure 2 shows the performance on up to 16 processors 
as obtained by comparing the wall-clock time with the single processor run. 


4 NUMERICAL EXAMPLES 

The methodology described above has been validated by comparison with 
both analytical solutions and experiments in [5]. The computations pre- 
sented below have been performed in order to study the feasibility of the 
method for large scale computations, its parallel performance, the eventual 
problems raised by such computations, as well as to gain some initial phys- 
ical insight into the influence of the fuselage and/or other surfaces on the 
sound field of a nacelle. Only inlet noise radiation of a single mode is com- 
puted in the cases presented here. The incoming pressure perturbation is 
specified as a circular duct eigensolution [19], 

p'{x, r, 6, t) = A mil J m (k mtt r) cos (k x x + mO + </> mM - ut) (15) 


where m is the azimuthal order of the spinning mode, J m the Bessel func- 
tion of the first kind, (j> mfl is the phase and A m/l the modal amplitude. The 
density and velocity perturbations are obtained from the equations of mo- 
tion. To keep the paper at a reasonable length, results are presented only 
for a stationary medium. Mean flow effects will be discussed in a subsequent 
work. 

Remarks on non-dimensionalization. In all cases, lengths were non- 
dimensionalized using the radius of the inlet duct as reference value. The 
same amplitude of the incoming acoustic wave was used for all computations. 
The pressure amplitude of the incoming acoustic modes at the fan face was 
set to A mii = 141.6 N/m 2 , and the sound pressure level (SPL) is reported in 
decibels relative to the usual reference pressure p re f = 2 x 10~ 5 N/m 2 . 


4.1 NACELLE ALONE 

The first computation was performed for the sound field radiated from an 
axisymmetric bell-mouth nacelle. The radius of the leading edge of the 
bell-mouth has been chosen to be 1/4 of the radius of the inlet duct, and 
no center-body (hub) was added. Radiation of the (2, 1) spinning mode 
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at reduced frequency oj = 8.0 in a quiescent fluid has been modeled on a 
grid made of E = 9919 elements. The polynomial degree used was N = 6, 
leading to an average of about five points per wavelength (largest element 
edge length is 1.06). There are a total of 3,402,217 Gauss-Legendre points 
in the domain. The nacelle casing extends on the x-axis (the symmetry axis) 
from x = —1.25 to x = 3.0, and the domain extends from x = —15 to x = 8, 
with a damping layer width of 3 in each direction. Incoming acoustic waves 
are specified on the source surface which is the circular disk with center at 
(0, 0, 0) and unitary radius. 

Figure 3 shows the SPL contours on the surface of the nacelle and in a 
plane that cuts the nacelle through its symmetry axis. Since no scattering 
surfaces other than the engine exist, the radiation is symmetric with respect 
to the engine axis. Figure 4 shows the contours in the plane z = —4. 
This computation was also used as an initial test of the radiation boundary 
conditions, which perform quite well in this case (very little waviness of the 
SPL contours is noticeable towards the interior limits of the damping layer). 
The main radiation lobe makes an angle of about 8 = 46 degrees with the 
x-axis, a value significantly different from 6 = 57 degrees obtained using the 
analytical formula derived by Rice [20] for infinitely thin cylindrical ducts, 


cos 9 = 



(16) 


where k mfl = 6.70613 is the circular duct eigenvalue for this mode. This 
difference is most probably caused by the large thickness of the nacelle lip 
(half of the inlet duct radius). 


4.2 FUSELAGE-NACELLE CONFIGURATION 

The second case considered a configuration obtained by joining the same 
nacelle with a cylindrical fuselage through an elliptic cylinder pylon. The 
cylinder representing the fuselage has a radius of 4.5 and its center is sit- 
uated at y = —7. The grid is made up of E = 17370 elements of the 
same order N = 6 (5,957,910 Gauss-Legendre points). The minimum and 
maximum values of the coordinates are (-18, -7, -9) and (8,10,8) respectively, 
and the damping layer has the same width as in the previous case, with 
the difference that the plane y = — 7 does not have an adjacent damping 
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layer in the direction perpendicular on it. Instead, a zero normal acoustic 
velocity boundary condition is used on this plane to simulate the effect of 
a symmetric configuration with engines in opposite phase on both parts of 
the fuselage. 

As can be seen in Figure 5, the presence of the fuselage creates a series 
of secondary radiation lobes due to the reflection of the nacelle lobes from 
the surface and to the various interferences that occur. The structure of the 
radiated sound field becomes thus quite complex. The presence and extent 
of the footprint of the SPL = 0 contour level in this plane is a measure 
of the effect of the fuselage in increasing the noise levels compared to the 
nacelle alone case. The figure also shows the trace of mesh decomposition 
on the solid surfaces defining the geometry. 


4.3 FUSELAGE- WING-NACELLE CONFIGURATION 

The geometry of the previous case was augmented with a wing of elliptic 
cross-section, slightly swept and mounted below the nacelle. The center 
of the ellipse defining the cross-section of the wing is situated in the plane 
z = —3. A grid of E = 22843 elements was generated in this case (for N = 6 
this corresponds to 7,835,149 discretization points). In all other respects, the 
definition of the domain is similar to the previous case. Shown in Figure 6 
are the SPL contours on the solid surfaces defining the geometry for this 
case. The presence of the wing is responsible for a complex, considerably 
asymmetric, sound field due to sound scattering from its leading and trailing 
edges. A clear picture of the effect of the wing can be obtained from figure 7, 
which shows the SPL levels in the plane y = 0. In this plane, the leading and 
trailing edge of the wing are situated at x = —7.2 and x = —1.5, respectively. 
As can be expected, although in the plane z = 0 (above the wing) the noise 
levels are, at most locations, higher than in the preceding case, in a plane 
below the wing the noise levels are remarkably lower, due to the upward 
reflection of the sound field by the wing surface, the difference immediately 
below the wing lying in the 10 — 15dB range. The corresponding contour 
plots are given in figures 8 and 9. It can be noticed in these figures that the 
damping layer region is not large enough to eliminate all reflections from 
the boundaries, as the contours are jagged towards the limits of the domain. 
The directivity at twelve nacelle radii from the center of the source disc for 
this case is compared with the nacelle alone data in figure 10. The shielding 
effect of the wing below the nacelle (negative 6 values) is clearly visible. 

For a quantitative comparison of the three cases, SPL data was extracted 
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in the field along several lines. Two such cases, representative for the in- 
creasing complexity of the sound field, are presented in figures 11 and 12. 
The interference patterns and various scattering effects that can be iden- 
tified from these plots render good predictions based on simple analytical 
methods extremely difficult, if not impossible. 


5 CONCLUSIONS 

The results presented here show that a more complete and accurate mod- 
eling of tone noise radiated from aircraft engines, including the effect of 
scattering from wing and fuselage, is now possible on the large distributed 
memory parallel computers available nowadays. This can have major im- 
pact in the design stage, in terms of affecting the noise footprint of the air- 
craft by changing the relative engine/wing positioning and altering the wing 
planform. The paper discusses a time domain methodology that is being 
developed to perform such studies. The reduced frequency for which results 
are presented here is lower than the usual values encountered at approach 
flight in current turbofans. However, it is our opinion that the paralleliza- 
tion strategy adopted for this method makes the study of the usual range 
of blade passage frequency values feasible on the supercomputers available 
nowadays. Further work, including the effect of mean flow and the sound 
radiated from the bypass exhaust, will be presented in the near future. 
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Figure 1: A ( = const plane through the master element, showing the 
position of the Gauss-Legendre points inside the element and the position 
of the flux points on the element faces. 
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Processors 



Figure 2: Actual wall-clock time parallel performance f< 
up to 16 processors: ideal (full line), 17370-elenrent grid 
884-element grid on SP3 (+). 







Figure 3: SPL contours for radiation of spinning mode (2, 1) at u = 8 from 
an axisymnretric bell-mouth nacelle alone. Nacelle surface and plane z = 0 
through its axis shown. Contour increment A = 2. 
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Figure 4: SPL contours for radiation of spinning mode (2, 1) at u) = 8 from 
an axisymmetric bell-mouth nacelle alone. Plane z = —4 shown. Trace of 
mesh decomposition shown on nacelle surface. 
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Figure 5: SPL contours for the radiation of spinning mode (2,1) at lo = 8 
in the fuselage-pylon-nacelle configuration in plane z = —4. Trace of mesh 
decomposition shown on the solid surfaces. Contour increment A = 4. 
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Figure 6: SPL contours for the radiation of spinning mode (2,1) at lo = 8 in 
the complete configuration with wing. Solid surfaces shown only. Contour 
increment A = 4. 
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Figure 7: SPL contours for the radiation of spinning mode (2,1) at u) = 8 in 
the complete configuration with wing. Plane y = 0 shown only. 
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Figure 8: SPL contours for the radiation of spinning mode (2,1) at lu = 8 
in the complete configuration with wing. Plane z = 0 and trace of mesh 
decomposition on solid surfaces shown. 
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Figure 9: SPL contours for the radiation of spinning mode (2,1) at lu = 8 
in the complete configuration with wing. Plane z = —4 and trace of mesh 
decomposition on solid surfaces shown. 
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Figure 10: Directivity in the plane y = 0 for the complete configuration with 
wing (dashed line) compared with the nacelle alone case (full line). 
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